Path-integral simulation of ice Ih: The effect of pressure 



(N 

o 

(N 



Carlos P. Herrero and Rafael Ramirez 
Instituto de Ciencia de Materiales de Madrid, Consejo Superior de Investigaciones 
Cientificas (CSIC), Campus de Cantoblanco, 28049 Madrid, Spain 
(Dated: January 5, 2012) 

The effect of pressure on structural and thermodynamic properties of ice Ih has been studied 
by means of path-integral molecular dynamics simulations at temperatures between 50 and 300 K. 
Interatomic interactions were modeled by using the effective q-TIP4P/F potential for flexible water. 
Positive (compression) and negative (tension) pressures have been considered, which allowed us to 
approach the limits for the mechanical stability of this solid water phase. We have studied the pres- 
sure dependence of the crystal volume, bulk modulus, interatomic distances, atomic derealization, 
and kinetic energy. The spinodal point at both negative and positive pressures is derived from the 
vanishing of the bulk modulus. For P < 0, the spinodal pressure changes from -1.38 to -0.73 GPa in 
the range from 50 to 300 K. At positive pressure the spinodal is associated to ice amorphization, and 
at low temperatures is found between 1.1 and 1.3 GPa. Quantum nuclear effects cause a reduction 
of the metastability region of ice Ih. 

PACS numbers: 31.70.Ks, 81.40.Vw, 65.40.De, 71.15.Pd 
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I. INTRODUCTION 

Despite the large volume of experimental and theoret- 
ical work on condensed phases of water, some of their 
properties still lack a complete understanding. This 
is mainly due to the peculiar structure of liquid and 
solid water, where hydrogen bonds between contiguous 
molecules give rise to properties somewhat different than 
those of most known liquids and solids (the so-called wa- 
ter "anomalies" ) fi _ — 

In the temperature region between 80 K and 250 K, ice 
Ih is the stable phase of water up to a pressure around 
0.2 GPa (that increases for increasing temperature). This 
range may seem rather narrow for some purposes, espe- 
cially if one compares it with other types of materials 
different from molecular solids. However, this pressure 
region is large if one takes into account the nature of the 
intermolecular interaction in water, which gives rise to 
the presence of hydrogen bonds. The region of mechan- 
ical stability of ice Ih, in which this phase is metastablc 
has been studied both experimentally and by computer 
simulations. At 80 K, simulations^ have predicted it to 
be metastable down to a negative pressure (tension) of 
about -1.4 GPa, whereas for positive pressure (compres- 
sion) its range of metastability is known to extend up 
to around 1 GPa, where ice Ih transforms into an amor- 
phous phase (the so-called high-density amorphous ice)£ 
This pressure-induced amorphization of ice Ih is an issue 
that has received in recent years a considerable amount 
of attention, along with the existence of different amor- 
phous phases which clearly differ in density.— ~— 

Computer simulation of condensed phases of water at 
an atomic level has a long history, dating back to around 
1970. Since then, important efforts have been focused 
on the development and refinement of empirical poten- 
tials to describe both liquid and solid phases of water, so 
that nowadays a large variety of this kind of potentials 
can be found in the literature J2r—. Many of them as- 



sume a rigid geometry for the water molecule, and some 
others include molecular flexibility either with harmonic 
or anharmonic OH stretches. In recent years, moreover, 
simulations of water using ab initio density functional 
theory (DFT) have been carried outJ£~— Nevertheless, 
it turns out that the hydrogen bonds in condensed phases 
of water are difficult to describe with presently available 
energy functionals, which causes that some properties 
cannot be accurately reproduced by DFT calculations. 21 
Some promising advances to improve the description of 
van der Waals interactions in water within the DFT for- 
malism have been recently presented %2u2L 

A limitation of ab-initio electronic-structure calcula- 
tions is that they usually treat atomic nuclei as classical 
particles, not including quantum effects like zero-point 
motion. These effects can be taken into account by us- 
ing harmonic or quasiharmonic approximations for the 
nuclear motion, but the precision of these approaches 
is not easily estimated when large anharmonicities are 
present, as can be the case for light atoms like hydrogen. 
To take into account the quantum character of atomic 
nuclei, the path-integral molecular dynamics (or Monte 
Carlo) approach has proved to be very useful, since in this 
procedure the nuclear degrees of freedom can be quan- 
tized in an efficient manner, thus including both quantum 
and thermal fluctuations in many-body systems at finite 
temperatures i 25 ' 26 Thus, a powerful method can consist 
in combining DFT to determine the electronic structure 
and path integrals to describe the quantum motion of 
atomic nucleii 18 ' 20 However, this procedure requires com- 
puter resources that would restrict enormously the num- 
ber of state points that can be considered in actual cal- 
culations. 

The phase diagram of water is now known up to tem- 
peratures and pressures on the order of 1000 K and hun- 
dreds of GPa, respectively. Present-day interatomic po- 
tentials are able to predict rather accurately the range of 
stability corresponding to the different known phases. 27 
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It is however important to estimate the influence of quan- 
tum nuclear effects on the stability range of the different 
ice phases. In this line, it is known that such quantum 
effects change the melting temperature T m of ice Ih at 
ambient pressure by some degrees, as manifested by the 
isotope effect on T m £& This kind of quantum effects has 
been less studied for liquid and solid water under exter- 
nal pressure. In particular, the spinodal line for ice Ih 
at negative pressures (which defines the limit of mechan- 
ical stability of this water phase) has been calculated by 
using classical molecular dynamics simulations. 5 Earlier 
studies of ice Ih using path- integral simulations have been 
carried out by using mainly effective potentials, and were 
focused on structural and dynamic properties of the solid 

Quantum nuclear effects become more relevant for light 
atomic masses, and are expected to be especially impor- 
tant in the case of hydrogen. Then, we pose the question 
of how quantum effects associated to the lightest atom 
can influence the structural properties of a solid water 
phase such as ice Ih, and in particular if these effects 
are relevant or detectable for the solid at different den- 
sities, i.e. under different external pressures. This refers 
to the crystal volume and interatomic distances, but also 
to the mechanical stability of the solid. In this context, 
it is usually assumed that increasing quantum fluctua- 
tions monotonically enhances the exploration of the en- 
ergy landscape. However, in certain regimes an increase 
in quantum fluctuations can lead to dynamical arrest, as 
shown for glass formation. 33 A similar mechanism could 
also retard the transition from crystalline to disordered 
(amorphous) phases. 

In the present paper we study ice Ih by path-integral 
molecular dynamics (PIMD) simulations at different 
pressures and temperatures, to assess the range of me- 
chanical stability of this water phase, as well as to an- 
alyze its structural properties. Interatomic interactions 
are described by the flexible q-TIP4P/F model, which 
has been recently developed and was employed to carry 
out PIMD simulations of liquid water 3 ^ and ice Ih at 
ambient pressure ! 28 ' 35 

The paper is organized as follows. In Sec. II, we de- 
scribe the computational method and the model em- 
ployed in our calculations. Our results are presented in 
Sec. Ill, dealing with the pressure dependence of crystal 
volume, bulk modulus, interatomic distances, atomic de- 
localization, and kinetic energy of ice Ih. Sec. IV includes 
a summary of the main results. 



II. COMPUTATIONAL METHOD 

We employ here the PIMD method to obtain equilib- 
rium properties of ice Ih at different temperatures and 
pressures. This method is based on an isomorphism be- 
tween the actual quantum system and a classical one, 
that appears after a discretization of the quantum density 
matrix along cyclic pathsi 36 ' 37 This isomorphism is in fact 



obtained by replacing each quantum particle by a ring 
polymer consisting of L (Trotter number) classical par- 
ticles, connected by harmonic springs with temperature- 
and mass-dependent force constant. Details on this sim- 
ulation technique can be found elsewherei 25 ' 26 We note 
that the dynamics employed in this method is fictitious 
and does not correspond to the actual quantum dynam- 
ics of the real particles under consideration, but it is 
used to effectively sample the many-body configuration 
space, giving precise results for the equilibrium proper- 
ties of the quantum system. An alternative way to derive 
equilibrium properties is the use of Monte Carlo sam- 
pling, but this procedure has been found to require for 
the present problem more computer resources than the 
PIMD method. A particular advantage of the latter is 
that in this case the codes can be more readily paral- 
lelized, an important factor for efficient use of modern 
computer architectures. 

Simulations of ice Ih have been carried out here in the 
isothermal-isobaric NPT ensemble (A, number of par- 
ticles; P, pressure; T, temperature), which allows us to 
obtain the equilibrium volume of the solid at given pres- 
sure and temperature. We have employed effective algo- 
rithms for performing PIMD simulations in this statis- 
tical ensemble, as those described in the literature^ 8 - - — 
Sampling of the configuration space has been carried out 
at temperatures between 50 K and 300 K, and pressures 
in the region of mechanical stability of ice Ih (which at 
50 K corresponds to the range from -1.4 to 1.1 GPa). 
Note that we have employed both negative (tension) and 
positive (compression) pressures in the simulations. For 
comparison with results of PIMD simulations, some sim- 
ulations of ice Ih were also carried out in the classical 
limit, which is obtained in our path integral procedure 
by setting the Trotter number L = 1. 

Our simulations were carried out on ice Ih supercclls 
with periodic boundary conditions. To check the in- 
fluence of finite-size effects on the results, we consid- 
ered orthorhombic supercells of two different sizes. The 
smaller one included 96 water molecules and had pa- 
rameters (3et, 2-^/30, 2c), where a and c are the stan- 
dard hexagonal lattice parameters of ice Ih, whereas the 
larger supercell included 288 molecules and had param- 
eters (4a, 3y/3a, 3c). Results obtained for both types of 
supercells coincided within the error bars due to the sta- 
tistical uncertainty associated to the simulation method. 
The flexibility of the simulation cell was taken into ac- 
count in the NPT simulations by treating the modules 
of the orthorhombic cell vectors as independent dynamic 
variables. In the considered supercells, and before the 
PIMD simulations, proton-disordered ice structures were 
generated by a Monte Carlo procedure, in such a way 
that each oxygen atom had two chemically bonded and 
two H-bonded hydrogen atoms, and with a cell dipole 
moment close to zero.— 

The interatomic interactions have been modeled by the 
point charge, flexible q-TIP4P/F model, recently devel- 
oped to study liquid water j2i and that was subsequently 



employed to study several properties of ice ; 28 ' 35 as well 
as water clusters 42 In line with the arguments discussed 
by Habershon et al.,— there are some reasons to take this 
model potential for the present study. First, most previ- 
ous works that have considered quantum effects in water 
phases have employed empirical potential models that 
were parameterized on the basis of earlier classical sim- 
ulations. Then, quantum simulations using those mod- 
els lead to "double counting" of quantum effects. Sec- 
ond, many of the empirical potentials previously used for 
quantum simulations of condensed phases of water treat 
H2O molecules as rigid bodiesj 30 ' 31 ' 43 This can be conve- 
nient for computational efficiency, but it disregards the 
important role of intramolecular flexibility in the struc- 
ture, dynamics, and thermodynamics of the condensed 
water phases^ Third, the significant anharmonicity of 
the O-H vibration of the water molecule is taken into 
account by anharmonic stretches in the q-TIP4P/F po- 
tential, vs. the harmonic potentials employed in the ma- 
jority of simulations that considered quantum effects. 

Technical details on the simulations presented here are 
the same as those employed and described in Refs. l28ll35l . 
In particular, the Trotter number L was taken propor- 
tional to the inverse temperature (L oc 1/T), so that LT 
= 6000 K, which allows one to keep roughly a constant 
precision in the PIMD results at different temperatures. 
Also, the time step At associated to the calculation of 
interatomic forces was taken in the range between 0.1 
and 0.3 fs, which was found to provide adequate conver- 
gence for the variables studied here. For given tempera- 
ture and pressure, a typical simulation run consisted of 
5 x 10 4 PIMD steps for system equilibration, followed 
by 6 x 10 5 steps for the calculation of ensemble average 
properties. At ambient pressure, the interatomic poten- 
tial q-TIP4P/F predicts melting of ice Ih at 251 but 
here we study this water phase up to 300 K, a temper- 
ature at which it was found to be metastable along our 
simulations. 

From PIMD simulations one can obtain insight into 
the atomic derealization at finite temperatures. This 
includes a thermal (classical) derealization, as well as 
a derealization associated to the quantum character of 
the atomic nuclei, which is quantified by the extension of 
the quantum paths in the path integral formalism. For 
each quantum path of a given particle, one can define the 
center-of-gravity (centroid) as 
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where r.; is the position of bead i in the associated ring 
polymer. Then, the mean-square displacement A 2 of the 
atomic nuclei (H or O) along a PIMD simulation run is 
defined as 
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FIG. 1: (Color online) Molar volume of ice Ih as a function of 
pressure, as derived from PIMD simulations at various tem- 
peratures: 50 K (squares), 150 K (circles), and 251 K (trian- 
gles). Error bars are in the order of the symbol size. Lines are 
guides to the eye. A solid diamond represents experimental 
data from Refs. I441I451 at atmospheric pressure and 273 K. 



where (...) indicates a thermal average at a given temper- 
ature. In connection with the quantum derealization of 
a particle, in this formalism it is interesting the spread of 
the paths associated to the particle, which can be mea- 
sured by the mean-square "radius-of-gyration" Q 2 of the 
ring polymers: 




(3) 



III. RESULTS 



A. Volume 



(2) 



We first present results for the equilibrium volume of 
ice Ih, as derived from our PIMD simulations at given 
pressure and temperature, in the pressure region where 
the solid turned out to be mechanically stable. These 
results are summarized in Fig. 1, where we show the 
pressure dependence of the molar volume at three dif- 
ferent temperatures: 50, 100, and 251 K. There are sev- 
eral important observations to make with regard to this 
figure. At a given temperature, we observe the usual vol- 
ume decrease for rising pressure, i.e., dV/dP < 0. In 
most of the pressure region considered, the slope dV/dP 
becomes less negative as the pressure is raised, which 
means d 2 V/dP 2 > 0. However, at about P = 0.5 GPa 
we observe a change in the trend of the first derivative 
as indicated by the presence of an inflection point with 
d 2 V/dP 2 — 0. At higher pressures, the slope increases in 
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absolute value until reaching the limit of mechanical sta- 
bility of the solid, close to P = 1 GPa. This change in the 
trend of the P — V curve, with the appearance of an in- 
flection point seems to be related to the proximity of the 
amorphization of the solid, with an important reduction 
in the molar volume£>&^ In fact, this amorphization is 
associated with a divergence in the compressibility of ice 
Ih, where the derivative dV/dP should diverge to — oo 
(see below). 

It is also remarkable that the volume-pressure curves 
cross at P ~ 0.2-0.3 GPa. At negative pressures we find 
that the larger volumes correspond to the higher tem- 
peratures, as expected for the usual thermal expansion 
of solids. However, at P > one finds just the opposite 
trend: the higher volume corresponds to the lower tem- 
perature, a fact associated to the negative thermal expan- 
sion occurring for ice Ih at low T. The potential TIP4P/F 
is able to reproduce this negative thermal expansion at 
atmospheric pressure, as discussed elsewhere^ At P = 
1 atm (10~ 4 GPa), the line corresponding to 251 K in 
Fig. 1 is still above those found for lower temperatures, 
but lies below them at higher (but relatively low) pos- 
itive pressures. This anomalous behavior of the crystal 
volume as a function of temperature is due to the nega- 
tive sign of the Griineisen parameter for TA vibrational 
modes42, For a mode in the n'th phonon branch with 
wave vector q, this parameter 7 n (q) is defined from the 
logarithmic derivative of its frequency with respect to the 
crystal volume 4£ 



Negative values of 7„(q) for TA modes in ice have 
been related to the tetrahedral coordination of water 
molecules, and a negative thermal expansion has been 
also observed in other solids with similar structures. In 
particular, it is well known for crystals with diamond and 
zinc blende structured In connection with this, nega- 
tive values of 7„ (q) for TA modes were also found to be 
related with the pressure-induced amorphization of ice 
IhiSlnfact, they cause that the relations between elastic 
constants necessary for crystal stability are violated at a 
certain applied pressure, giving rise to mechanical insta- 
bility of the solid. In this line, instabilities in transverse 
acoustical modes in ice Ih were also found at low tem- 
peratures from incoherent inelastic neutron scattering^ 
For P < 0, the absolute value of the slope dV/dP in- 
creases fast, and eventually diverges for a finite value of 
the pressure, which depends on temperature. This is as- 
sociated with a divergence of the compressibility of the 
solid, which gives the limit for the mechanical stability of 
this ice phase (spinodal point) at negative pressures. On 
the other side, for positive pressures the stability limit at 
P ~ 1 GPa is related to the amorphization of the solid. 
In fact, at these pressures we find in the PIMD simula- 
tions a sudden reduction in the volume of the solid, ac- 
companied by the breaking and reordering of hydrogen 
bonds (see below). Note that the possibility of studying a 



solid phase in metastable conditions, close to a spinodal 
line is limited by the appearance of nucleation events, 
which cause the breakdown of the crystal structure and 
the surge of the equilibrium phase. In this kind of atom- 
istic simulations the probability of such nucleation events 
at low temperatures is relatively low, and the metastable 
range of the solid that can be explored is rather large. In 
fact, one can go to conditions close to the spinodal lines 
at positive and negative pressures, especially at P < 0. 

At 250 K ice Ih was found to be stable along our simu- 
lations up to P = 0.9 GPa. At such pressures, the H-bond 
network of ice Ih collapses and the crystal transforms into 
the amorphous phase with an important volume decrease. 
At T = 251 K and P = 1 GPa, this decrease amounts 
to about 18% of the crystal volume. Something similar 
happens at lower temperatures, and at 100 K we found 
a volume reduction of 24% in the amorphization process. 
This is in line with the pressure-induced amorphization 
of ice Ih first observed at T = 77 K, and occurring at 
about P = 1 GPa£ This water phase is the so-called 
high-density amorphous ice.— — Note that at higher T 
(in the order of 250 K) amorphization of ice Ih has not 
been experimentally observed, and a transition to ice III 
at P ~ 0.2 GPa could be expected. Such transitions 
between crystalline ice phases are usually not observed 
from the kind of atomistic simulations employed here, 
and other types of simulations employing thermodynamic 
integration would be required for this purposed The im- 
portant point here is that we obtain at each temperature 
the pressure range in which ice Ih is metastable with the 
employed potential model, i.e., we are not discussing the 
transitions between equilibrium (crystalline) phases, as 
they would appear in the phase diagram of water. 

Our results are in line with those derived from clas- 
sical molecular dynamics simulations. 5,51 In particular, 
Sciortino et al. 5 employed the TIP4P potential, and 
found that the solid suffers reversible changes at 77 K 
for compression (P > 0) up to ~ 1.5 GPa. At higher 
pressures, ice was found to become amorphous in an ir- 
reversible way, i.e. it remained in the amorphous phase 
when decreasing the pressure to low values in the order 
of 1 atm. 

The molar volume of ice Ih has been determined exper- 
imentally in different ways, e.g., x-ray, optical, mechan- 
ical, calorimetric, acoustical, or nuclear methods (see 
Ref. [H and references therein). Measurements of dif- 
ferent authors typically deviate from each other by up to 
about 0.3%d Our results at atmospheric pressure are in 
line with those of Refs. fTlTuii considered in Ref. HH the 
most accurate determinations at normal pressure and T 
— 273 K. These data are represented in Fig. 1 by a solid 
diamond. 



B. Bulk modulus and spinodal pressure 

The compressibility of liquid water and ice shows pecu- 
liar properties associated to the hydrogen-bond network 



present in these condensed phases. It turns out that their 
compressibility is smaller than what one could expect at 
first sight from the large cavities present in their struc- 
ture, which could probably collapse under pressure with- 
out water molecules approaching close enough to repel 
each other. However, this is not the case, and for ice 
Ih in particular the hydrogen bonds holding the crystal 
structure are known to be rather stable, as revealed by 
the relatively high pressure necessary to break down the 
H-bond network^ 

The isothermal compressibility n of ice, or its in- 
verse the bulk modulus [B = 1/k = — V(dP/dV)r] 
can be straightforwardly obtained from our PIMD sim- 
ulations in the isothermal-isobaric ensemble. In this 
ensemble the isothermal bulk modulus can be calcu- 
lated from the mean-square fluctuations of the volume, 
a v — (^ 2 ) — (^) 2 i by employing the expressio: 



,53.54 
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k B T(V) 
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ks being Boltzmann's constant. This expression was 
used earlier to obtain the bulk modulus of different types 
of solids from path-integral simulations ! 1 

At atmospheric pressure we have found results for the 
bulk modulus of ice Ih similar to those derived from ex- 
periments. In particular, at T = 251 and 300 K, we found 
B = 9.1 GPa and 8.4 GPa, respectively. There appears 
in the literature a dispersion of data for the isothermal 
compressibility (or bulk modulus) of ice Ih at the melting 
temperature and normal pressure. The values of B ob- 
tained here are in line with B = 8.49 GPa at T = 273 K, 
corresponding to the equation of state derived by Feister 
and Wagner— for this water phase from a set of exper- 
imental data. It is difficult to estimate the precision of 
this value of the bulk modulus due to the lack of error 
bars for the published data, and the uncertainty may be 
large, as suggested by the dispersion in experimental data 
obtained by different authors (see Ref. [52| and references 
therein). The uncertainty is even larger at lower temper- 
atures, where we could not find any direct experimental 
data for the isothermal compressibility. The bulk modu- 
lus calculated from our PIMD simulations decreases for 
rising temperature (a fact commonly found in solids) , and 
this behavior appears to be true for all studied pressures, 
except for P > 1 GPa, where this general trend seems to 
be inverted (see below). 

We will first discuss our results for negative pressure 
(solid under tension). For P < 0, the bulk modulus 
decreases as the pressure becomes more negative (see 
Fig. 2), which means dB/dP > at all temperatures 
under consideration. Then, B is expected to vanish at a 
temperature-dependent pressure P s , which corresponds 
to the spinodal point defining the limit of mechanical 
stability of ice Ih at negative pressures (where the com- 
pressibility diverges to infinity). 

At a given temperature and near the spinodal pressure, 
our results for the bulk modulus can be fitted to the ex- 
pression B ~ (P — Pg) 1 / 2 , which allows us to obtain P s 
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FIG. 2: (Color online) Pressure dependence of the bulk mod- 
ulus of ice Ih as obtained from quantum PIMD simulations at 
several temperatures: 50 K (squares), 150 K (circles), and 251 
K (triangles). Error bars show the statistical uncertainty in 
the values of B found from the simulations. Lines are guides 
to the eye. 



from a linear fit of B 2 vs P. This pressure dependence of 
B can be understood in the following way^ At a tem- 
perature T, and calling F the free energy, the spinodal 
point is given by the condition d 2 F/dV 2 \v B — 0, which 

0. Thus, close to P s one 



can be rewritten as dP/dV\v 3 
can write along an isotherm: 

P = P S + c 2 (v - V s f + C 3 (V - v s f + 



(6) 



with constants C'2 and C3 independent of V. Assuming 
C2 / (which seems to be true here), one has B « 
2y/U^V s {P - Ps) 1 / 2 for small values of the difference P - 
P s (for discussions on the case C2 = 0, see Ref. [57l - l59| ). 
Such a pressure dependence for B, or its equivalent: 



P-Ps 



{V~Vsf 



(7) 



has been obtained earlier for ice and Si02 cristobalite, as 
well as for rare-gas solids close to their stability limitsi^i^ 
This means in our case that the exponents describing the 
singularity are consistent with those of the mean-field 
spinodal^ given by Eq. ([7]), where P s and V s are the 
values of the pressure and volume on the spinodal. Thus, 
we have obtained P s at each temperature by fitting our 
results for the bulk modulus at P < to an expression 
B ~ (P — Pg) 1 / 2 . The results for P s , as derived from our 
PIMD simulations, are shown in Fig. 3 by solid circles. 
P s becomes less negative as the temperature is raised, 
going from -1.38 GPa at 50 K to -0.73 GPa at 300 K. 

We now turn to the results at positive pressures (com- 
pression). For P > 0, the bulk modulus increases mod- 
erately until reaching a maximum value at about 0.3-0.5 
GPa. This maximum value of B is lower for higher T, 
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FIG. 3: (Color online) Calculated spinodal pressure for ice Ih 
as a function of temperature. Solid circles are data obtained 
from PIMD simulations, whereas open squares represent data 
yielded by classical simulations. Labels 'q' and 'cl' indicate 
quantum and classical, respectively. A solid diamond cor- 
responds to the experimental point for amorphization of ice 
Ih given in Ref. @ (denoted here as MI84). Labels indicate 
the regions in the P — T diagram where one finds the differ- 
ent phases: ice Ih, high-density amorphous (HDA), and gas. 
Lines are guides to the eye. 



as shown in Fig. 2, and appears to be slightly shifted to 
lower pressures as temperature increases. Strassle et al&> 
measured the equation of state of ice Ih at 145 K under 
hydrostatic pressure up to 0.5 GPa, and fitted their re- 
sults for the pressure dependence of the volume to a Mur- 
naghan equation. They obtained at zero pressure a bulk 
modulus B — 9.85(47) GPa, somewhat lower than that 
derived from our simulations at 150 K: B = 11.9 GPa. 

At pressures larger than 0.5 GPa the bulk modulus de- 
rived from our PIMD simulations decreases rather fast, 
until reaching at each temperature a certain pressure at 
which B approaches zero and ice Ih becomes mechani- 
cally unstable. We will call P' s this spinodal pressure. 
Close to P' s one has a relation between B and P similar 
to that given above, but valid here for P' s > P, so that : 
B ~ (P' s — P) 1 / 2 . Fitting our results to this expression 
we find that the bulk modulus extrapolates to zero for P s 
= 1.11 and 1.25 GPa, at 50 and 100 K, respectively. At 
higher temperatures, this extrapolation cannot be done 
from our PIMD results in a reliable way, since the solid 
collapses into the amorphous phase at pressures clearly 
lower than that corresponding to the metastability limit 
of the crystal, where the bulk modulus is still far from 
vanishing. This is due to a kinetic phenomenon such 
as the nucleation of a new phase, which prevents close 
approach to the spinodal, as indicated above. Thus, at 
250 K ice Ih was metastable in our PIMD simulations 
only up to P = 0.9 GPa, where the bulk modulus is still 



larger than 7 GPa, but the ice crystal transforms into the 
amorphous phase with an appreciable volume decrease. 
Arguments equivalent to those presented here for the ice 
instability associated to a vanishing bulk modulus, have 
been given in connection with the elastic moduli of this 
hexagonal crystal^ 

In Fig. 3 we present values of the spinodal pressure P' a , 
as derived from our PIMD simulations, along with those 
corresponding to negative pressures, P s (represented as 
solid circles). In both cases the spinodal pressure in- 
creases as the temperature is raised. At positive pres- 
sure, we show only results up to 100 K, as at higher T 
the extrapolation employed to determine P' s is unreliable, 
as explained above. These results are compared with 
those derived from classical molecular dynamics simula- 
tions (open squares). At positive pressure, P' s derived 
from classical simulations is larger than that found from 
PIMD, and the opposite occurs for the spinodal P s at 
negative pressures. In both cases, quantum and clas- 
sical results are found to lie closer as temperature is 
raised. Altogether we find that quantum effects reduce 
the metastability range of ice Ih. This means that at 
low temperature the pressure region in which this water 
phase is metastable is reduced by about 0.3 GPa, both 
at negative and positive pressures. For comparison with 
the calculated spinodal, we also show in Fig. 3 the value 
of P = 1 GPa given by Mishima et alr± for the laboratory 
amorphization pressure at 77 K (solid diamond). Spin- 
odal pressures for ice Ih have been calculated earlier by 
Sciortino et al^ from classical molecular dynamics sim- 
ulations with the TIP4P potential. These authors found 
at low temperatures values close to the results of our 
classical simulations, but at high T they obtained values 
somewhat higher (less negative) than those found here, 
e.g., at 250 K they found P s — —0.75 GPa vs our value 
of -0.99 GPa. This difference is presumably due to dif- 
ferences between the interatomic potentials employed in 
both calculations, as the TIP4P deals with rigid water 
molecules and the q-TIP4P/F takes into account their 
flexibility, allowing for bending and stretching of the in- 
tramolecular bonds. Simulations of the amorphization of 
ice Ih using classical molecular dynamics were also car- 
ried out by Tse and Klein ; 63 ' 64 who found ice amorphiza- 
tion at pressures around 1.2—1.3 GPa at temperatures 
between 80 and 100 K, close to the spinodal P' s obtained 
here. 

We remark that values of the bulk modulus calculated 
from PIMD simulations in the isothermal-isobaric ensem- 
ble show relative error bars larger than those obtained for 
other variables (e.g., molar volume, kinetic energy, or in- 
teratomic distances), as a consequence of the statistical 
uncertainty in the volume fluctuations try, employed to 
calculate the bulk modulus. An alternative way to de- 
rive B consists in calculating numerically the derivative 
dV/dP from the P—V equation of state at temperature 
T, We have checked that this method gives results that 
agree within error bars with those derived from the vol- 
ume fluctuations through Eq. However, the calcula- 
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FIG. 4: (Color online) Mean distance between oxygen atoms 
in nearest-neighbor water molecules in ice Ih, as a function of 
pressure. Symbols represent results of PIMD simulations at 
different temperatures: 50 K (squares), 150 K (circles), and 
251 K (triangles). Lines are guides to the eye. The dashed line 
represents the results of classical simulations at 50 K. Solid 
circles represent O-O distances derived from the structural 
data obtained by Strassle et aZ.— from neutron diffraction 
experiments of D2O ice Ih at 145 K. 

tion based on the pressure derivative of the volume may 
be less reliable in the pressure region where B changes 
fast, as happens close to P = 1 GPa (see Fig. 2). 

C. Interatomic distances 

In this section we show results for interatomic dis- 
tances in ice Ih between atoms in the same and adjacent 
molecules. This can shed light on the structural changes 
suffered by the crystal when temperature and/or pres- 
sure are modified. We first show in Fig. 4 the mean dis- 
tance between oxygen atoms in neighboring molecules as 
a function of external pressure. Open symbols represent 
results of PIMD simulations at different temperatures. In 
each case, we give values of d(O-O) in the region where 
the crystal was found to be metastable in the PIMD sim- 
ulations. This includes, as in previous figures, positive 
and negative values of P. For increasing pressure, d(0- 
O) is reduced, as could be expected from the associated 
volume decrease (dV/dP < 0). Nevertheless, contrary to 
the case of the crystal volume, the isothermal lines cor- 
responding to the 0-0 distance do not cross each other, 
so that at higher temperature one finds a larger distance 
in the whole metastability region of ice Ih. 

For an homogeneous and isotropic expansion (or con- 
traction) of the crystal, one expects for small volume 
changes: AV/V = 3A//Z, I being any distance in the 
solid, and one would expect this relation to be fulfilled, in 



particular, for the intermolecular distances in the solid. 
However, although this relation is roughly followed in 
the parameter range considered here, it is obvious that 
it cannot be strictly accomplished, given the apparent 
different trends of the curves shown for crystal volume 
and 0-0 distance in Figs. 1 and 4. In particular, vol- 
ume curves at different temperatures cross each other, 
whereas those corresponding to the 0-0 distance follow 
a regular pattern without crossings. The molar volume 
at low temperature (T = 50 K) decreases by 7.0% when 
going from atmospheric pressure to P = 1 GPa, from 
where one would expect a relative change in the average 
0-0 distance of 2.3%. This value is somewhat higher 
than the change actually obtained from the simulations, 
i.e., 1.8%. On the other side, at negative pressure we 
find a relative volume increase of 20.4% at P = -1.35 
GPa (close to the spinodal pressure P s ), which would 
correspond to a relative change in the 0-0 distance of 
6.8%. In fact, we find from the simulations an increase in 
the 0-0 distance of 6.3%, slightly smaller than that ex- 
pected for an homogeneous and isotropic crystal expan- 
sion. This indicates that expanding or contracting the 
ice crystal is accompanied by a distortion of the tetrahe- 
dral network of water molecules. This kind of distortion 
is also related to the negative thermal expansion of ice 
Ih occurring at low temperatures, and that was earlier 
found to be reproduced by PIMD simulations with the 
q-TIP4P/F potential^ 

For comparison with the results of PIMD simulations, 
we also show in Fig. 4 results of classical molecular dy- 
namics simulations with the same interatomic potential, 
which are shown for T = 50 K by a dashed line. At at- 
mospheric pressure the average 0-0 distance obtained in 
the quantum simulations is 0.032 A larger than that de- 
rived from classical simulations, which means a difference 
of 1.2%. This difference increases at negative pressures, 
and at P — — 1.2 GPa and T = 50 K we find an increase 
of 0.069 A due to quantum effects, i.e., a 2.4% of the 
0-0 distance. These numbers are consistent with the 
volume increase associated to quantum effects. In fact, 
at P = 1 atm the quantum result for the molar volume is 
3.4% larger than the classical one, vs an increase of 7.2% 
obtained at P = —1.2 GPa (see Ref. l35l for results at am- 
bient pressure and several temperatures). In Fig. 4 we 
also display 0-0 distances derived from the structural 
data obtained by Strassle et alr±^ from neutron diffrac- 
tion experiments of D2O ice Ih at 145 K (solid circles). 
Note that 0-0 distances in ice may change with the hy- 
drogen isotope^ so that these data for D2O ice at 0.02 
and 0.25 GPa could differ from those corresponding to 
H2O ice, but in any case they are useful for comparison 
with our PIMD results. 

We now consider the interatomic distances inside wa- 
ter molecules in the ice crystal, as derived from our sim- 
ulations at different pressures. In Fig. 5 we have plot- 
ted the mean intramolecular O-H distance as a function 
of pressure at three temperatures. Concerning this fig- 
ure, there appear two results that should be emphasized. 
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FIG. 5: (Color online) Pressure dependence of the mean in- 
tramolecular O-H distance, as derived from simulations of ice 
Ih. Open symbols represent results of PIMD simulations at 
different temperatures: 50 K (squares), 150 K (circles), and 
251 K (triangles). Error bars are in the order of the symbol 
size. Lines are guides to the eye. 
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FIG. 6: (Color online) Pressure dependence of the spatial de- 
localization of hydrogen nuclei (protons) in ice Ih, as derived 
from PIMD simulations. Symbols indicate the mean-square 
displacement at several temperatures: 50 K (squares), 100 
K (circles), 150 K (triangles), and 251 K (diamonds). Lines 
are guides to the eye. The dashed line is the mean-square 
displacement of hydrogen derived from classical molecular dy- 
namics simulations at 50 K. 



First, at a given pressure we observe that the O-H dis- 
tance decreases as temperature is raised, contrary to the 
usual expansion of atomic bonds for increasing T. Sec- 
ond, at a given temperature, we find that the O-H bond 
distance increases for rising pressure, which may seem 
also surprising given the pressure-induced volume con- 
traction. Both facts are due to the characteristic struc- 
ture of ice with hydrogen bonds connecting contiguous 
water molecules, and can be explained by the same ar- 
gument: An increase in the intramolecular O-H distance 
is associated to a weakening of the O-H bond, and is 
caused by a hardening of the intermolecular H bond. For 
rising T, molecular motion is enhanced, so that hydro- 
gen bonds become softer, and the average intermolecular 
O-H distance increases (in line with the rise in 0-0 dis- 
tance shown above and in Ref. Ha ). This gives rise to 
a strengthening of the intramolecular O-H bonds, with 
the associated decrease in interatomic distance in water 
molecules. This hardening of intramolecular O-H bonds 
for rising temperature has been observed experimentally 
and reported in the literature^ For increasing pressure 
at a given temperature, we find something similar: the 
intermolecular 0-0 distance is reduced as the pressure 
rises (see Fig. 4), in such a way that the intermolecu- 
lar H bonds become stronger, causing weaker and longer 
intramolecular O-H bonds. The rise in O-H distance de- 
rived from our simulations in the range between and 
1 GPa is consistent with an increase of about 2 x 10~ 3 
A/GPa, given in Refs. [66ll67l for high-pressure ice phases, 
although smaller values have been found in Ref. l68l . 

Note, however, that the change of intermolecular 0-0 



distance with pressure is much larger than that of the 
intramolecular O-H distance. Thus, at 50 K the change 
in d(0-0) in the whole metastability region of ice Ih 
(from P = -1.38 to 1.12 GPa) is about 8%, whereas the 
change in intramolecular O-H distance is about 1%. This 
is indeed a consequence of the direct effect of pressure 
in changing the intermolecular distance, to be compared 
with an indirect effect on the intramolecular O-H bond 
strength caused by intermolecular H bonds. 

We finally note that the intramolecular O-H distance 
increases appreciably due to quantum nuclear effects. In 
fact, a classical simulation at T — 50 K and ambient 
pressure yields a mean O-H distance of 0.969 A, to be 
compared with 0.984 A derived from PIMD simulations, 
which means an increase of 1.5% in the bond length. This 
difference is much larger than the temperature-induced 
change in d(O-H) shown in Fig. 5. However, the increase 
due to quantum motion is rather constant in the whole 
pressure range studied here, i.e. at negative and positive 
pressures we find a rise in O-H distance of about 1.5%. 

D. Atomic delocalization 

Here we present results for the atomic mean-square dis- 
placements defined in Eq. One expects that this 
spatial delocalization will be larger for hydrogen than 
for oxygen, due to the smaller mass of the former. In 
Fig. 6 we display values of Aj; for hydrogen in ice Ih as a 
function of pressure. Data were derived from PIMD sim- 
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FIG. 7: (Color online) Pressure dependence of the spatial 
derealization of oxygen nuclei, as derived from PIMD simu- 
lations. Symbols represent the mean-square displacement 
at several temperatures: 50 K (squares), 100 K (circles), 150 
K (triangles), and 251 K (diamonds). Lines are guides to the 
eye. The dashed line is the mean-square displacement of oxy- 
gen derived from classical molecular dynamics simulations at 
50 K. 



ulations at different temperatures. At a given pressure, 
A 2 increases as the temperature is raised, in the whole 
region where the solid is found to be metastable. At 50 
K (the lowest temperature shown here), A 2 changes very 
little as a function of pressure, except close to P = 1 GPa, 
where it increases by a small amount as one approaches 
the spinodal point (which at this temperature was found 
to be P' s = 1.12 GPa; see above). At higher tempera- 
tures, changes in A 2 as a function of pressure are more 
appreciable. At all temperatures considered, we find an 
increase in the atomic derealization when approaching 
one of the spinodal points, at positive or negative pres- 
sure. The increase at positive pressure (close to amor- 
phization of the solid) is clearly larger than that found 
at negative pressures. In Fig. 6 we also present the mean- 
square displacement of hydrogen atoms in Ih, as derived 
from classical simulations at 50 K (dashed line). These 
results are clearly lower than those of the PIMD simula- 
tions at the same temperature. Thus, at ambient pres- 
sure the classical result is A 2 = 2.68 x 10 -2 A 2 , to be 
compared with A 2 = 7.88 x 10 -2 A 2 found in the quan- 
tum simulations. 

For oxygen we find mean-square displacements that 
change with pressure in a way similar to those of hydro- 
gen. This is shown in Fig. 7 at the same temperatures 
as in Fig. 6 for hydrogen. At atmospheric pressure, we 
find for oxygen A 2 = 4.08 x 10~ 2 A 2 at 50 K vs 0.175 
A 2 at 251 K. Values of A 2 are larger for hydrogen than 
for oxygen, but the difference between both decreases for 
increasing temperature. Thus, at T — 50 K and P 



1 atm, A 2 for hydrogen is 1.93 times the mean-square 
displacement for oxygen, whereas this ratio decreases to 
1.13 at 251 K. This is in line with a larger quantum delo- 
calization for hydrogen, caused by its lower mass, and the 
effect of the mass becomes less relevant as temperature 
increases. 

Note that the mean-square displacement A 2 of a given 
atomic nucleus can be decomposed into a contribution 
Q 2 coming from the spread of the quantum paths [see 
Eq. and another one, C 2 , given by the spatial dis- 
placement of the centroid r of the paths i i 35 At 50 K 
and 1 atm, Q 2 represents a 62% of A 2 for hydrogen and 
a 25% in the case of oxygen, which reflects the fact that 
the quantum contribution to the atomic derealization is 
more important in the case of hydrogen. This can also be 
seen by directly comparing values of Q 2 for the different 
atomic species, which result to be 4.8 times larger for H 
than for O at 50 K and atmospheric pressure. 

At low temperatures, the quantum contribution 
(spreading of the paths) , Q 2 , dominates the spatial der- 
ealization of hydrogen, since the centroid displacement, 
C 2 , converges to zero as T — > K. The opposite happens 
in the high-temperature limit (unreachable here for sta- 
bility reasons), where the quantum contribution Q 2 even- 
tually disappears, as corresponds to the classical limit. 
We note in passing that at low temperatures the quantum 
paths associated to hydrogen have an average extension 
of about 0.15 A, much smaller than the H-H distance in 
a water molecule, thus justifying the neglect of quantum 
exchange between protons in the PIMD simulations. 

In Fig. 7 we also show mean-square displacements of 
oxygen atoms, as derived from classical simulations at 
50 K (dashed line). As one could expect, these values 
are smaller than those derived from PIMD simulations, 
but in the case of oxygen the difference between both 
sets of results, although nonnegligible, is not so impor- 
tant as for hydrogen. In fact, at ambient pressure we 
found for oxygen A 2 = 2.75 x 10 -2 A 2 in classical sim- 
ulations vs 4.07 x 10~ 2 A 2 derived from PIMD simula- 
tions. Note that the classical result is similar to that 
found for hydrogen at the same pressure and tempera- 
ture, but the quantum A 2 for H is clearly larger than 
that of O atoms. Also, quantum effects on the mean- 
square displacement for oxygen decrease rapidly as tem- 
perature is raised. Thus, at T — 150 K and P = 1 atm 
we find in the classical simulations A 2 = 0.091 A 2 , to 
be compared with A 2 = 0.102 A 2 derived from PIMD 
simulations, i.e., quantum effects increase A 2 by a 12% 
at this temperature. 



E. Kinetic energy 

The kinetic energy of atomic nuclei in a solid or 
molecule depends on the mass and derealization of the 
considered nucleus. Thus, although in a classical ap- 
proach, each degree of freedom contributes to the kinetic 
energy by an amount that depends only on temperature 
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FIG. 8: (Color online) Kinetic energy of hydrogen in ice Ih 
as a function of pressure. Open symbols indicate results de- 
rived from PIMD simulations at different temperatures: 50 K 
(squares), 150 K (circles), and 251 K (triangles). Error bars 
are less than the symbol size. Lines are guides to the eye. 



FIG. 9: (Color online) Kinetic energy of oxygen atoms in ice 
Ih as a function of pressure. Open symbols indicate results 
derived from PIMD simulations at different temperatures: 50 
K (squares), 100 K (circles), 150 K (triangles), and 251 K 
(diamonds). Error bars are less than the symbol size. Lines 
are guides to the eye. 



(feflT/2, as given by the equipartition principle), in a 
quantum approach the kinetic energy, Ef., gives infor- 
mation on the environment and interatomic interactions 
seen by the considered particle. In particular, a typical 
quantum effect related to the atomic motion in solids is 
that the kinetic energy at low temperature converges to 
a finite value associated to zero-point motion, contrary 
to the classical result where Ek vanishes at K. Path 
integral simulations allow us to calculate the kinetic en- 
ergy of the quantum particles under consideration, which 
is related to the spread of the quantum paths. In fact, 
for a particle with a certain mass and at a given tem- 
perature, the larger the mean-square radius-of-gyration 
of the paths, Q^, the smaller the kinetic energy, in line 
with the expectancy that a larger quantum derealiza- 
tion causes a reduction in the kinetic energyj 25 ' 35 Here we 
have calculated Ek by using the so-called virial estimator, 
which has an associated statistical uncertainty apprecia- 
bly lower than the potential energy of the system i 39 i 69 

We present separately the kinetic energy of oxygen and 
hydrogen atoms in ice Ih. In Fig. 8 we display Ek for 
hydrogen as a function of pressure at three different tem- 
peratures, as derived from our PIMD simulations. At 
each considered temperature, Ek increases as pressure 
rises, corresponding to an overall increase of vibrational 
frequencies. Ek also increases with temperature, but the 
data for T — 50 and 100 K (not shown in the figure) 
are almost indistinguishable, since at these temperatures 
the hydrogen vibrations are nearly in their ground state. 
Higher vibrational modes are excited at higher tempera- 
tures, and at 251 K we observe an increase in Ek larger 
than 0.1 kJ/mol in the whole metastability range of ice 



Ih. Values of Ek calculated at P — 1 atm coincide with 
those presented earlier from PIMD simulations in the 
range from 210 to 290 

The kinetic energy of hydrogen atoms in ice Ih was ob- 
tained by Rcitcr et alJ^ at ambient pressure from inelas- 
tic neutron scattering experiments. These authors found 
at 269 K a kinetic energy Ek — 12.7 kJ/mol, somewhat 
lower than those found here at P — 1 atm in the tem- 
perature range from 250 to 290 K, i.e., Ek for hydrogen 
between 14.39 and 14.44 kJ/mol. 

In Fig. 9 we show the pressure dependence of the ki- 
netic energy of oxygen in ice Ih at four temperatures. 
As for hydrogen, for oxygen Ek also grows for increas- 
ing pressure at a given temperature. However, now the 
isothermal lines at 50 and 100 K are clearly distinct by 
more than 0.2 kJ/mol, as a consequence of the higher 
mass of oxygen, which produces a lower vibrational en- 
ergy and an appreciable excitation of vibrational modes 
at 100 K. From the results for the kinetic energy of H and 
O shown in Figs. 8 and 9, we observe that at atmospheric 
pressure and low temperatures, Ek for a hydrogen atom 
in ice Ih is about 4 times larger than that of an oxygen 
atom. At 50 K the kinetic energy of hydrogen increases 
from 13.96 kJ/mol at P = -1.35 GPa to 14.32 kJ/mol 
for 1 GPa. This represents an increase of 0.36 kJ/mol, 
i.e., 2.6% in the whole pressure range. For oxygen, the 
increase in Ek in the same region is about twice that of 
hydrogen (AEk — 0.75 kJ/mol), but its relative change 
is much more important for O (a 26.5%). 

Close to the amorphization pressure of ice we observe 
a slow increase in the kinetic energy for both hydrogen 
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and oxygen, without any appreciable change in its trend. 
This contrasts with the trend observed for the atomic 
mean-square displacements Aji shown above in Figs. 6 
and 7. This is a consequence of the fact that Ek is re- 
lated to the spread of the quantum paths, as measured 
by the mean-square displacements Q%, and not to the 
centroid motion, which represents displacements of the 
overall paths, irrespective of their size and shape. The 
fast increase in close to the spinodal points is mainly 
due to an increase in the centroid derealization 

IV. SUMMARY 

We have presented results of PIMD simulations of ice 
Ih in the isothermal-isobaric ensemble at different pres- 
sures and temperatures. This kind of simulations have 
allowed us to study this water phase in the pressure re- 
gion where it is metastable, and approach the spinodal 
points at which it becomes mechanically unstable. At 
P < we found a spinodal pressure P 8 ranging from -1.4 
to -0.7 GPa in the temperature range from 50 to 300 K. 
At positive pressures, we obtained a spinodal in the range 
P' s ~ 1.1-1.3 GPa, that is well defined at temperatures 
T < 100 K. At higher temperatures, kinetic processes fa- 
vor amorphization of the solid for pressures clearly lower 
than the estimated spinodal pressure. 

Although the molar volume decreases with increasing 
pressure, the interatomic distances in the solid change 
in a peculiar way that remembers their temperature de- 
pendence. Thus, the distance between oxygen atoms in 
neighboring molecules decreases for increasing pressure, 
but the intramolecular O-H distance becomes larger for 
higher pressure. This is a consequence of the hydrogen 
bonds connecting contiguous molecules, which become 
stronger as the volume (or 0-0 distance) is reduced, 
causing a weakening of the intramolecular O-H bonds. 

For the bulk modulus of ice Ih, we obtain a maximum 
at 0.3-0.5 GPa, which increases slowly as the tempera- 
ture is lowered. At higher and lower (negative) pressures, 



B is found to decrease and extrapolates to zero at the 
corresponding spinodal pressure (P s or P s '). Close to the 
spinodals, we observe an increase in the atomic dereal- 
ization for both oxygen and hydrogen atoms, especially 
at temperatures in the order of the melting temperature 
of ice Ih. For the kinetic energy, however, we do not ob- 
serve any anomalous effect, and it is found to increase 
smoothly as temperature or pressure is raised. 

We have assessed the importance of quantum effects 
by comparing results obtained from PIMD simulations 
with those yielded by classical simulations. Concerning 
the stability of ice Ih, we have found that quantum ef- 
fects reduce the metastability region of this solid phase, 
for both positive and negative pressures. At low tem- 
peratures, this means that both spinodal pressures P s 
and P' s are shifted by about 0.3 GPa. Structural vari- 
ables also change when quantum nuclear motion is con- 
sidered. Thus, the crystal volume, interatomic distances, 
and atomic derealization suffer appreciable modifica- 
tions in the range of temperature and pressure considered 
here. Most important is the increase in the mean-square 
displacement of hydrogen, which contributes to weaken 
the intermolecular H bonds, and to strengthen the in- 
tramolecular O-H bonds. 

An extension of this work can consist in studying the 
high-density amorphous phase obtained from ice Ih un- 
der pressure, as well as other amorphous water phases, 
which may show peculiar properties as a function of pres- 
sure and temperature. Path-integral molecular dynamics 
simulations with interatomic potentials similar to that 
employed here can be also adequate to study these solid 
phases. 
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